/*  scfcsqrt.cpp  by K.Tsuru  */
// function ID = 9101
/**************************************************
SComplex class
It provides the sqrt(z) = x + i*y.
It is assumed that -pi <= arg(z) <= pi, then x >= 0
and y has the same sign as the imaginary part of z.
I referred to gcc's "complext.cc".
--------------------------------------------------
z = a +i * b
= (r, theta) ... Representation in polar coodinate.
sqrt(z) = (sqrt(r), theta/2) = x+i*y

Let theta = t.
a = r cos(t), b = r sint.
cos(t/2) = sqrt{ [1+cos(t)] / 2} =sqrt{[1 + (a/r)] / 2}
then we get
x = sqrt(r)cos(t/2) = sqrt{(r+a)/2} > 0 (always positive).
In the same way
y = [+|-]sqrt(r)sin(t/2) = [+|-]sqrt{(r-a)/2}

(x+iy)^2=x^2-y^2+2ixy=a+ib
x^2-y^2=a, 2xy=b
x>0 then y has the same sign of b.
**************************************************/
#ifndef SN_H
#include "sn.h"
#endif
SComplex Csqrt(const SComplex& z){
    if(z.Imag().Sign() == 0) { // b = 0, contains the case of a = 0  ver. 2.18
        SDouble r = Sqrt( Dabs( z.Real() ) );
        if(z.Real().Sign() > 0) return SComplex(r, 0.0);
        return SComplex(0.0, r);
    }
    
    SDouble r = Cabs(z), x, y, a = z.Real(), b = z.Imag(); // z = a + i b
    if(a.Sign() > 0) { // a > 0, -pi<theta<pi
        SDouble rx = RecSqrt( DsDiv(r + a, 2) );    // rx = sqrt{2/(r+a)} = 1/x
        x = DReciprocal(rx);                        // x = 1/rx = sqrt{(r+a)/2} > 0
        y = DsDiv(b, 2) * rx;                       // y = (b/2)*(1/x) = b/(2*x) the same sign of b
    } else {    // a <= 0, r - a > 0, pi<theta<2*pi
        SDouble ry = RecSqrt( DsDiv(r - a, 2) );    // ry = sqrt{2/(r-a)} = 1/y
        y = DReciprocal(ry);                        // y = 1/ry = sqrt{(r-a)/2}
        x = DsDiv(b, 2) * ry;                       // x = (b/2)*(1/y) = b/(2*y)
        if(b.Sign() < 0){
            x.ChangeSign();  y.ChangeSign();        // x = -x = |x| y = -y <0 --> xy=b/2 < 0
        }
    }
    return SComplex(x, y);
}